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SUMMARY 


This report describes the details of the work performed under Element A, "Im- 
proved Numerical Methods for Turbulent Viscous Recirculating Flows," of NASA 
sponsored Hot Section Technology (HOST) program. The objective in this study 
is to develop accurate and efficient numerical schemes for predicting complex 
flows . 

Five discretization schemes were selected for preliminary evaluation based on 
simple flows. These schemes were evaluated on the basis of accuracy, stabil- 
ity, and ease of extension to multidimensions. This initial evaluation led to 
the selection of three schemes — linear flux-spline, cubic flux-spline, and con- 
trolled numerical diffusion with internal feedback (CONDIF) for further evalua- 
tion in two-dimensional flows. The accuracy was assessed by solving a series 
of test problems that included scalar transport, laminar flows, and turbulent 
flows. The numerical results were compared with analytical solutions, experi- 
mental data, and fine grid numerical solutions. For these test problems, the 
linear flux-spline scheme was superior to other schemes and was selected for 
incorporation into a computer program for three-dimensional flows. This scheme 
was further evaluated by solving three-dimensional flows. 

To improve the computational efficiency, a coupled solution approach for the 
fluid flow equations was adopted. This algorithm was incorporated in conjunc- 
tion with the linear flux-spline scheme. For two-dimensional flows, this ap- 
proach led to a factor of 2-3 reduction in execution times compared with the 
sequential algorithms. To extend the coupled solution approach to three— dimen- 
sions, a plane-by-plane solution strategy was followed. For sample flows, such 
a procedure was robust and fast convergent. 
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I . INTRODUCTION 


1.1 OBJECTIVES 

The overall objective of the hot section technology (HOST) aerothermal modeling 
program — Phase II is to improve the accuracy of the current aerothermal models 
for gas turbine combustors. Specifically Element A, Improved Numerical Methods 
for Turbulent Viscous Recirculating Flows, seeks improvements in the accuracy 
of the differencing schemes and the computational efficiency of the solution 
algorithms. Improvements in these areas would allow accurate predictions of 
complex flows in a cost-effective manner without requiring an excessively large 
number of grid points. Further, the use of improved differencing schemes, 
which do not suffer significantly from the false diffusion problem, will allow 
an objective evaluation of various models of gas turbine combustion processes. 

1.2 APPROACH 

The approach followed in this study was to select a number of discretization 
schemes and solution algorithms and assess them on the basis of accuracy, com- 
putational efficiency, stability, and ease of extension to multidimensions. 

For this preliminary evaluation, test problems with known analytical solutions 
were selected. Based on the results of these test problems, three discretiza- 
tion schemes were chosen for an in-depth evaluation for two-dimensional flows. 
In addition, a coupled solution approach for the fluid flow equations was con- 
sidered. 

The selected differencing schemes were incorporated into computer programs for 
two-dimensional flows. The accuracy characteristics of these schemes were as- 
sessed by solving a series of test problems, which included scalar transport, 
laminar flows, and turbulent flows. The scheme with superior performance was 
then combined with the coupled solution approach for the fluid flow equations. 
For some sample flows, this coupled solution approach was compared with an 
iterative sequential algorithm. 

The selected discretization scheme and solution algorithm were then incorpor- 
ated into a computer program for three-dimensional flows. The improvements in 
accuracy and computational efficiency were demonstrated by solving some sample 
three-dimensional test problems. 

1.3 OUTLINE OF THE REPORT 

In Chapter II, an overview of the finite-volume procedure is presented. This 
is followed by a brief description of the selected differencing schemes and 
the coupled solution strategy. The results of these schemes for sample test 
problems are presented. The schemes are evaluated on the basis of various cri- 
teria mentioned earlier. 

In Chapter III, the three differencing schemes with superior performance are 
described in detail. In the development, the expressions have been given for 
one coordinate direction only; expressions in other directions can easily be 
obtained by simple axis transformation. 

In Chapter IV, the results of two- and three-dimensional test cases are pre- 
sented. For two-dimensional problems, the numerical results have been compared 
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with available analytical, experimental, and fine-grid numerical solutions. 

The scheme with superior performance in two-dimensional situations was used to 
solve three-dimensional problems. 

In Chapter V, a coupled solution approach for fluid flow equations is de- 
scribed. For sample flows the performance of such an approach has been com- 
pared with that of an iterative procedure. 

In Chapter VI, the present work is briefly reviewed and suggestions for further 
work are given. 
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II. DESCRIPTION OF THE SELECTED DISCRETIZATION SCHEMES 
AND SOLUTION ALGORITHM 


2.1 FINITE VOLUME METHOD 

This section includes preliminary details of the finite- (control-) volume 
method. The main objective of this discussion is to prepare a common frame- 
work within which various schemes, to be discussed later, can be accommodated. 
For ease of presentation, the development has been restricted to two-dimen- 
sional steady flows in Cartesian coordinates. 

The conservation equation for a dependent variable $ can b e expressed in the 
following general form (Ref 1*): 

3Jx 3JY _ o 
3x + 9y 

where 

Jx = P U * " r ^ 

jy = pv<i> - r |* 

T is the "effective” diffusion coefficient and S is the source of 4>. 

Integrating equation (1) over the control volume around the grid point P (see 
Figure 1) gives: 

J e “ J w + J n ~ J x = (Sc + S P 4>p) 

where the source term S has been linearized. The quantities J e , J w , J n , and J s 
are integrated fluxes over the control-volume faces. That is, 

J e = /Jx dy over the interface e 

= <pu«, - r jg> e 4y <5) 

To estimate the value of the flux at the control-volume faces, approximations 
for the convected value and the gradient (8$/dx) e are needed. These two 
approximations are the essence of any differencing scheme. A differencing 
scheme, via a profile assumption for $, leads to an expression for the flux 
at a control-volume interface that involves the values of 4> at grid points 
in the vicinity of the face. When the expressions for fluxes are substituted 
in equation (4), an algebraic equation involving $'s is obtained. This equa- 
tion can be expressed in a general form as: 

ap 4»p = Z a^ $nb + bfi + b s (6) 


( 1 ) 

( 2 ) 

(3) 


♦References are listed at the end of this report. 
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Figure 1. Control volume for the two-dimensional situation. 


where nb refers to the neighbors of the grid point P under consideration (nb = 
E, W, N, S), b s includes the physical sources, and bfi is an apparent source. 
The expression for b f ^ involves either distant neighbors along a coordinate 
line or adjacent lines or fluxes at control volume interfaces. 

2.2 DISCRETIZATION SCHEMES 

The discussion of the discretization schemes selected in this study is mostly 
qualitative. The mathematical details have been kept to a minimum. Instead, 
emphasis has been placed on the features of each scheme that result in better 
accuracy, the computational molecule involved, and the factors leading to spa- 
tial oscillations and the possibility of lack of convergence. 

The following five discretization schemes were selected for a preliminary as- 
sessment of their accuracy: 

1. locally analytic differencing scheme (LOADS) 

2. linear flux-spline scheme 

3. cubic flux-spline scheme 

4. a modified central differencing scheme, CONDIF 

5. mass weighted skew upwind differencing scheme (MWSUDS) 

The Power-law differencing scheme (PLDS) has been chosen as the representative 
of all lower-order schemes based on the one-dimensional convection-diffusion 
equation without a source. 

2.2.1 Power-Law Differencing Scheme (PLDS) 

The Power-law differencing scheme (Ref 1) is based on a curve fit to the exact 
solution of one-dimensional convection-diffusion equation in the absence of a 
source. For the situation shown in Figure 2, the final discretization equation 
may be expressed as: 
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ap $p = a E 4> E + a w 4> w + b (7) 

where 

a E = D e A (|P e |) + max [-F e , 0] (8a) 

a W = &w A (|P W |) + max [F w , 0] (8b) 

ap = a E + ayj — Sp Ax (8c) 

b = S c Ax (8d) 

F = (pU) (8e) 

D = T/6x (8f ) 

and 

P = F/D (8g) 

The function A( P ) is defined as: 

A ( IP I ) = max [0, (1 - 0.1IPI) 5 ] (9) 


where max [ ] stands for the largest of the quantities contained within it. 

Note that by a proper choice of A (1PI), other schemes can also be recovered. 
For the commonly used schemes, function A ( |P|) is listed in Table I. 


Table I. 


The function A ( P) 

Scheme 

Upwind 

Hybrid 

Exponential 

Power-law 

Central 


for various schemes. 


A ( P ) 

1 

max [ 0 , 1-0 . 5 1 P I ] 

P /[exp(IPI) -1] 
max [ 0, <1-0.1 IP I ) 5 ] 
1-0.5 IPI 


U 



Figure 2. Grid configuration for a one-dimensional situation. 
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Since the Power-law differencing scheme is based on a purely one-dimensional 
flux balance, it is expected to perform well in the regions in which the flow 
is aligned along the grid lines and in which convection is balanced primarily 
by streamwise diffusion rather than cross-stream diffusion or sources. If such 
idealized conditions are not encountered, the locally one -dimensional assump- 
tion gives rise to significant numerical errors (false diffusion). The im- 
proved schemes attempt to include the effect of flow-to-grid skewness, lateral 
transport, and the source in the ^-distribution between the grid points. 

2.2.2 Locally Analytic Differencing Scheme (LOADS) 

LOADS (Ref 2) takes into account the influence of the lateral transport pro- 
cesses and source terms in the derivation of the ^-distribution between two 
grid nodes. 

The scheme is based on the solution of the equations: 

fc <p«+-r |J> - is - <pv*-r !*» = kx <io) 

* - r “> - i s - k ^ - r - * < u > 

where Kx and Ky can be viewed as source terms for flows in x and y directions, 
respectively. They take into account the actual source and the apparent source 
in a given direction due to the net transport in the other direction. The 
quantities Kx and Ky are evaluated from the numerical solution currently avail- 
able and are assumed to be constant between two grid points. This leads to a 
^-distribution along a direction (say x) of the form: 

$ = A + B exp (jTTc) + Cx (12) 

The coefficients in the final discretization equation are identical to those 
resulting from the exponential scheme except for an additional source that is 
a function of Peclet numbers in the x and y directions. 

The computational molecule for LOADS involves 5 points in two dimensions (7 
points in 3-D). The influence coefficients are always positive and the alge- 
braic equations are diagonally dominant. The source terms Kx and Ky, however, 
are calculated from nonconverged (incorrect) solutions and are treated explic- 
itly. Such a practice may lead to convergence difficulties, especially if the 
equations are strongly coupled. The presence of the apparent source may also 
result in spatial oscillations in the numerical solution. 

2.2.3 Linear Flux-Spline Scheme 

The lower-order schemes, such as hybrid and exponential schemes, represent the 
solution of the one-dimensional convection diffusion equation with zero source. 
These schemes, therefore, are based on the assumption of constant total (con- 
vection and diffusion) flux within a control volume and do not respond to the 
presence of sources or sinks and multidimensional effects. The linear flux- 
spline scheme (Ref 3) removes this shortcoming by asstiming a linear variation 
for the total flux, as shown in Figure 3. 
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Thus, the ^-profile in the x-direction for this scheme is the solution of 
the equation: 


d<t> 

- r s - J i + 


J ui J i 

Ax 


( 13 ) 


This can be written as: 

4> = A + B exp x) + Cx (14) 

where the coefficients A, B, and C can be expressed in terms of 4»i » J i» 
and J • 

The influence coefficients for the linear flux-spline scheme are identical to 
those from the exponential scheme except that the linear flux-variation leads 
to an additional source term. This source term involves the difference of 
fluxes at the control-volume faces and thus allows the scheme to respond to 
the presence of sources and multidimensional effects. 

The principles underlying the linear flux-spline and LOADS are very similar. 
They differ only in the manner in which the multidimensional effects are intro- 
duced in the ^-profile; i.e., the definition of the apparent source. In 
LOADS the contribution of transport processes in other directions is calcu- 
lated, within the iterative process, from the nonconverged flow and ^-fields. 
The linear flux-spline scheme seeks a field solution for the fluxes also. So, 
during the iterative process, the fluxes are related to each other and satisfy 
the conservation principle; therefore, the linear flux-spline scheme is ex- 
pected to be superior to LOADS. As with LOADS, the possibility of wiggles due 
to apparent source also exists in the linear flux-spline scheme. In addition, 
the details regarding the computational molecule are identical to those for 
LOADS . 
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2,2,4 Cubic Flux-Spline Scheme 


The cubic flux-spline scheme (Ref 4) is an improved variant of the linear flux- 
spline scheme. In this scheme, a cubic profile is assumed for the variation 
of the total flux of $ and the mass flux (pU) within a control volume. 

These assumptions lead to a ^“profile of the form: 

pH 2 

♦ = A + B exp ( i x) + Cx + Dx + Ex 3 (15) 


As for the linear flux-spline scheme, the cubic flux-spline scheme also in- 
volves additional "flux" source terms that allow this scheme to respond to the 
presence of the source and multidimensional effects. 

The rather sophisticated assumptions involved in the cubic flux-spline scheme 
make it computationally more complex than other schemes considered. The coef- 
ficients are no longer simple algebraic function of the Peclet number but in- 
volve numerical integration. In addition, field solution is required for the 
fluxes and their derivatives as well as the derivatives of the mass fluxes. 
Similar to the linear flux-spline scheme, this scheme also leads to a five- 
point formulation in two dimensions. The presence of apparent source may lead 
to wiggles. 

Although the cubic flux-spline scheme involves more computational effort than 
other schemes, this disadvantage may be offset by its higher order of accuracy. 
Thus for comparable accuracy, the cubic flux-spline scheme would require fewer 
of grid points. 


1*2 15 Controlled Numerical Diffusion with Internal Feedback (CONDin 

CONDIF (Ref 5) is a modified central dif ferencinge scheme (CDS) that retains 
the second-order accuracy of central differencing scheme but eliminates the 
overshoots and undershoots that occur if the grid Peclet number exceeds 2. 

The CONDIF modifies CDS by introducing a controlled amount of numerical diffu- 
sion based on local gradients. 


CONDIF uses central differencing for points at which the grid Peclet number is 
less than 2. For points where the Peclet number exceeds 2 and the dependent 
variable exhibits a monotonic behavior, the central differencing scheme is mod- 
ified as follows (for 1-D flows): 


*p«rVE*U ,b 

where 


t 

*E 


*E, CDS + 
V CDS + 


[ I (pu) e l + (pU) e ] 
5 

M(pU) w l- (pU) v ] 


+ Aw/R 
+ A e R 


A e = [ l(pU) e | + (pU) e ] / 4 


(16) 

(17a) 

(17b) 

(17c) 
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(17d) 


A w = t I(PU)„| - (pU) w ] / 4 
R = 04>/3x) e /(34>/9x) w 

In these expressions, CDS and a W, CDS are t ^ ie coefficients resulting 
from the central differencing scheme. At the nodes where the Peclet number 
exceeds 2 and the ^-profile is nonmonotonic, the upwind scheme is used. 

Since in most practical flows the dependent variables go through an extremum 
only at relatively few points, CONDIF uses the central differencing scheme in 
its original or modified form over most of the grid points. By using the CDS, 
the numerical diffusion is reduced. 


The computational molecule for CONDIF involves only 5 points in two dimensions. 
The coefficients in the modified CDS become nonlinear since they involve the 
gradients of the dependent variable. As a result, the coefficients for a 
linear problem need to be recomputed within the iterative process. However, 
the possibility of wiggles has been eliminated by adding controlled amount of 
numerical diffusion in the regions of steep gradient. Unlike the previous 
schemes, CONDIF does not involve any apparent source in the discretization 
equation. 


— Mass-Weighted Skew Upwind Differencing Scheme (MWSUDS) 

This scheme (Ref 6) seeks to improve the original skew upwind differencing 
scheme (SUDS) (Ref 7), by ensuring the positivity of the coefficients and thus 
eliminating the spatial oscillations. 

Before presenting the details of MWSUDS, the cause of negative coefficients in 
the skew upwind differencing scheme is analyzed. A typical control volume is 
shown in Figure 4. For illustration purposes, the variation of the south coef- 
ficient with the flow angle is examined. The value of $ s will appear in 
the convective fluxes through the east and south faces. The convective flux 
through the east face is given by: 

m e t e = m e t** 6 $S + (*• ~ Fxe) $p] ( 18 ) 

Similarly, the convective flux through the south face is: 

“s+s = m s 4>SW + (1* - Fys) 4> s ] ( 19 ) 

In above expressions, m e and m 8 are the mass flow rates through the 

east and south faces, respectively. Fxe and Fys are weighting factors defined 

as: 


Fxe 


V 6x/2 

r ~&r 


Fys 


U 6y/2 
s 

V 6x 
s 


( 20 ) 

( 21 ) 
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Figure A. A typical control volume arrangement. 


From equations (18) and (19), the convective part of the south coefficient a s 
can be obtained as: 

a s > -m e Fxe + m 8 (1 - Fys) (22) 

The coefficient becomes negative when 
m e Fxe > m s (1 - Fys) 


(23) 


This may lead to a lack of boundedness. Note that the presence of negative 
coefficient does not necessarily lead to wiggles in the numerical solution. 

In MWSUDS , the coefficients are forced to remain positive (nonnegative in 
strict sense) under all possible conditions. This is accomplished as follows: 

1. The weighting of $ s from the east face and south faces are li “ k * d ‘ 

This coupling is affected via the choice of velocity component V e needed 
for calculating the skew angle at the east face. In original SUDS, Vg 
is taken as the average of V at four neighboring points. In the modified 
skew scheme V e is approximated as V s . 

2. The weighting factors are expressed in terms of face-skewed mass flow rate 
components as : 


Fxe 


|pVgl Sx/2 

lpu e | 6y 


(24) 
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To ensure stability, the following condition should be satisfied: 

m e Fxe m s (1 - Fys) (25) 

Combining equations (24) and (25), the following relation can be obtained: 


0,5 i 1 - Fys 

(26) 

or Fys £ 0.5 

Therefore, an upper limit of 0.5 is placed on the skewness factor. 

The computational molecule for MWSUDS involves 9 points in two dimensions and 
27 points in three dimensions. Because of this, the tridiagonal matrix algor- 
ithm (TDMA) may not be a suitable solver for the algebraic equations. 

The modified scheme, like its original version, accounts for the flow-to-mesh 
skewness but not for the effect of sources and cross-stream diffusion. It re- 
duces to the upwind scheme for one-dimensional flows. 


2.3 PRELIMINARY EVALUATION OF THE SCHEMES 


The five selected schemes were used to solve a number of model one- and two- 
dimensional problems. In this section, the relative performance of these 
schemes is evaluated. The test problems considered are the following: 


1. uniform flow in a pipe with a heat source 

2. solid body rotation with logarithmic temperature distribution 

3. recirculating flow with a prescribed heat source 

In addition to the results from the selected schemes, results from two other 
popular schemes, QUICK (Ref 8) and SUDS (Ref 7) have also been presented for a 
relative assessment. These results have been taken from Reference 4. 


2.3.1 Uniform Flow in a Pipe with a Heat Source 

This problem consists of a uniform velocity flow in a pipe. A heat source 

S - 1 + 2x + II cos (IIx) + [II 2 sin (IIx)-2]/P (27) 

is prescribed such that the temperature distribution given by 

T ex = 1 + x + x 2 + exp (-P) exp (Px) + sin (Fix); 0 £ x £ 1 (28) 

represents an exact solution to the governing equation. 

31 _ 1 s!t = S (29) 

dx P 3x 2 

This problem was solved for a range of Peclet numbers (P) using 10 uniformly 
spaced control volumes in the computational domain. The results are presented 
in Figure 5 in terms of the percent average error, E av , defined by 

E = i I|(l- T /I )l* 100 (30) 

av N c ex' 
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where T c is the computed solution and N is the number of internal nodes. 

In the low Peclet number range all schemes result in comparable errors with 
the exception of the cubic flux-spline scheme, which leads to errors that are 
about two orders of magnitude smaller. Also, the errors for the linear flux- 
spline scheme, LOADS, and CONDIF are almost constant over the Peclet number 
range considered. Other schemes, however, show an increase in the error as 
the Peclet number exceeds 2. 

2.3.2 Solid Body Rotation with Logarithmic T emperature Distribution 

This test problem, proposed by Runchal (Ref 9), involves heat conduction in a 
fluid in solid body rotation. The geometry flow field for the problem is shown 
in Figure 6. This flow situation is one-dimensional in cylindrical coordi- 
nates. However, when solved in Cartesian coordinates, the problem becomes 
fully two-dimensional with finite convective and diffusive contributions. 

The governing transport equation for this problem is 


8x 


(P°T - 1 S 



(pvr - r . 


0 


(31) 


where U and V are the velocity components in the x and y directions, respec 
tively. For a solid body rotation, these are given by 

U = 2y, V = -2x 


For a uniform diffusivity, T, given by 

r = i/p, 


(33) 



Figure 6. The solid body rotation problem: the geometry and control volume. 
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the exact solution to the governing equations is obtained as: 


Tex 


2 2 
, In (x + v ) 

2 In 3 


( 34 ) 


This problem was solved using 10 x 10 uniformly spaced control volumes. The 
average error, E av , in the numerical solution for the various schemes con- 
sidered is shown in Figure 7. 

All improved schemes, with the exception of MWSDS, perform considerably better 
than the lower-order power law differencing scheme at high Peclet numbers. 

The cubic flux— spline scheme performs exceptionally well, the error being al- 
most two orders of magnitude smaller than that from other schemes. The perfor- 
mances of linear flux-spline, LOADS, CONDIF, and QUICK schemes are comparable 
over the Peclet number range considered. The use of MWSDS leads to results 
that are inferior to those from the Power— law scheme because of the inconsis- 
tency in the flux expressions at the control-volume faces. That is, the flux 
expression obtained for the east side of a control volume is not the same as 
the flux expression for the west side of the next control volume. The validity 
of the upwind mass flow approximation is also questionable for an arbitrary 
flow field. 


2_t3 , 3 Recirculating Flow with a Prescribed Heat Source 


In this test problem, transport of a scalar in a prescribed flow field is con- 
sidered. The velocity field represents a recirculating flow and is given by 


U = -Uy (1 - y2) (l _ x 2)2 


(35) 


V = +Ux (1 - x 2 ) (1 - y 2 ) 2 


(36) 


A source S, defined as 


S = P[4 - 2(x 2 + y 2 )] 


(37) 


is prescribed such that the exact temperature distribution is given by 

T ex = (1-x 2 ) (1-y 2 ), -1 i x, y i 1 ( 38 ) 

This problem was solved using a uniform grid with 15 control volumes in each 
direction. 


Figure 8 shows the maximum normalized error in the computed solution from vari- 
ous discretization schemes. For this problem, the flux-spline schemes give 
results with minimum numerical error. As with the previous problems, the cubic 
flux-spline is slightly more accurate than the linear flux-spline scheme. The 
performance of LOADS is not very different from that of the Power-law scheme 
(PLDS) . 

Convergence problems were encountered with LOADS for P 2. 100. This indicates 
the dominance of the apparent source, which was approximated in a rather simple 
manner . 
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Average error 







The results from CONDIF are not much better than those from the Power— law 
scheme, especially at high Peclet numbers. The use of MWSUDS results in lower 
numerical error compared with PLDS for P > 10. 

2.4 SUMMARY 

The selected convection-diffusion schemes were used to solve three test prob- 
lems for which analytical solutions are available. Based on the results of 
these problems, the following general conclusions can be drawn: 

1. The improved discretization schemes, with the exception of LOADS and 
MWSUDS, produce more accurate results compared with the PLDS for all 
problems considered. 

2. The simple treatment of the flux-terms in LOADS may result in lack of 
convergence, especially at higher Peclet numbers. The convergence behav- 
ior is expected to deteriorate as the coupling between various equations 
becomes stronger. 

3« For the test problems considered in this study, the performance of 
MWSUDS was not significantly superior to that of PLDS. As stated 
earlier, this behavior may be a consequence of the inconsistency in the 
discretized flux expressions at the control-volume faces. 

4. The performances of LOADS and MWSUDS appear to be problem dependent. 

5. Due to mild ^-distributions in the selected problems, none of the 
schemes produced spatial oscillations. 

On the basis of these findings and other studies, e.g., Ref 3, 4, 5, the fol- 
lowing three schemes were selected for an in-depth evaluation in two-dimen- 
sional flows: 

1. linear flux-spline scheme 

2. cubic flux-spline scheme 

3. CONDIF 

2.5 SOLUTION OF THE FLUID FLOW EQUATIONS 

The solution of the fluid flow equations in the primitive variable formulation 
requires special consideration due to the special role of pressure. The pres- 
sure appears in the momentum equations, but there is no pressure term in the 
continuity equation. Thus there is no equation for an explicit evaluation of 
the pressure field. The commonly used iterative methods based on SIMPLE (Ref 
1) or its variants rewrite the continuity equation as an equation for pressure 
or pressure correction and the equations are solved sequentially in a decoupled 
manner. The performance of such an approach depends on a proper choice of the 
underrelaxation factors and deteriorates as the number of grid points is in- 
creased. 

An alternative to the sequential SIMPLE-based approach is to obtain a direct 
solution to the whole set of continuity and momentum equation using a sparse 
matrix variant of Gaussian elimination. Such a method implicitly retains the 
coupling between various equations and eliminates the need for an equation for 
pressure (correction). This coupled approach has been used by Vanka and co- 
workers (Ref 10, 11) and Braaten (Ref 12) for solving various sample flows in 
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conjunction with a lower-order discretization scheme (hybrid and Power-law). 

In all these studies the Yale sparse matrix package (YSMP) (Ref 13) was used 
for the LU decomposition. The use of the direct solution strategy reduced the 
computational times by factors of 5—10 compared with the SIMPLE algorithm. In 
view of its impressive performance, the YSMP-based direct solution approach 
was selected for improving the computational efficiency of the overall solution 
algorithm. 
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III. DESCRIPTION OF SELECTED DISCRETIZATION SCHEMES 

In this section, the three selected convection-diffusion schemes are described 
in detail. Since the schemes are based on the solution of one -dimensional con- 
vection-diffusion equation, most details have been restricted to a single di- 
mension. The expressions for other directions can be obtained by an axis 
transformation. The derivations have been presented for a general nonuniform 
grid. By a proper specification of various geometrical parameters, the expres- 
sions are valid for a scalar as well as the velocity components. For fluid 
flow calculations, a staggered grid arrangement is used. 

3.1 LINEAR FLUX-SPLINE SCHEME 

3.1.1 One-Dimensional Convection-Diffusion 

The steady-state form of the one-dimensional convection-diffusion equation is 


dJ _ „ (39) 

dx " S 

where J is the total flux (convection and diffusion) of the dependent variable 
4>, and S is the source or sink of <j>. The total flux is given by 


J 


pu<t> - r 


d$ 

dx 


(AO) 


To obtain the variation of $ within a control volume, a profile assumption 
for J is required. The lower-order methods, such as the exponential scheme, 
assume that the total flux is uniform within a control volume. Such an assump- 
tion leads to a linear profile for <fr in the conduct ion- like problems and an 
exponential profile in convection-diffusion problems. Due to this assumption, 
the resulting schemes do not respond to the multidimensional effects and the 
presence of sources. The linear flux-spline scheme is based on the assumption 
that the total flux varies in a piecewise linear manner as shown in Figure 9. 
Thus for a control volume around the grid point i, the total flux is given by 

J = Ji + (Ji+i - Ji) (x/Axj) < 41 > 


J 1 ♦ 2 
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If the expression for J is substituted in equation (41), a differential equa- 
tion governing the variation of $ is obtained. The solution of the resulting 
equation requires the variation of the diffusivity I and the mass flux (pU). 
These quantities are assumed to be uniform over a control volume. With this 
assumption, the following expression for 4> results: 

♦i = A + Bx + C exp (pUx/D (42 

where 


A = 


B = 


C = 


(pu), 


+ ( A tr J i ) 

+ vpU)^ Ax^ 


(PU), 


(pU) i AXj 


(pU). 


11±1 
(pU) 


+ <^7) 


exp (-P^ 


P i " 


(pU) i Ax 1 


(43a) 

(43b) 

(43c) 

(43d) 


Equation (42) gives the variation of 4> within a control volume. For two ad- 
jacent control volumes, the 4>-profiles are such that they imply the same to- 
tal flux at the common interface. In addition, they must also give a unique 

value of 4> at the common interface. With reference to Figure 10, this condi- 
tion is: 


♦ 


+ 

i-1 


= 4* 


i 


(44) 


where and are defined as the values of 4> at the left and 

right interfaces, respectively, of the control volume around the grid point i. 

From equations (42) and (43) 


v. (Ax') 2 1 

1 J. - L — 


~ p i 1 AXj 

*i = *i 6 + A(-P i ) J i ' Ax, r7 (J i' J i+i ) G ( -v 


i ‘i 


(45) 


J 1 - 1 » 


♦l-l 4>1 

. -o — f o 

I 4 , I 


! J i + 1 


1 - i 


1 


1 


TE84-1680A 


Figure 10. Two adjacent control volumes. 
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. + . p l-l 1 to l-l T 

^i-1 “ M-l e “ + T J i 

A <- p iIi> i-i 




where 


( 46 ) 


o+ (P u >i-i *»j_i 


i-1 


(47) 


i-1 


A(P) = 


e P -l 


(48) 


G(P) = 


e (P-1) t 1 


(49) 


Combining equations (43) to (49), the expression for flux Jj can be written as 
J i = (Di - Ej ♦i) + B i (J r J i+1 ) + C i (^-J^) (50) 


where 


D i 

= H i exp (P^) 

(51) 

E i 

= Hj exp ( 

-V 

(52) 

B i 

(Ax") 2 

1 

- H. G (-P ) 
i 

(53) 

= Ax i r 

c i 

to i-X 

r ul B i G < p i-i> 

(54) 

and 




H i 

/to' 1 

+ to i-i 1 \ 

(55) 

II 

H* | 

> 

1 

P I> ’ r i-i *< p i!i>/ 


In equation (50), the expression (Dj - E^ <^j) is identical to that ob- 

tained from the lower-order exponential scheme. The terms involving Bj and Cj 
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are additional contributions of the flux-spline scheme based on the linear 
variation of the total flux. 

For ease of presentation, let 

Jli = Bj (Jj - Ji + i) + (Jj - Ji_i) (56) 

so that equation (50) may be expressed as: 

Ji = (Di $i-l - Ei 4>i) + Jli (57) 

In this derivation, the mass flux (pU) was assumed to be uniform within a 
control volume. To allow for variable mass flux, a correction term is added 
in the flux expression. The choice of this term is such that the governing 
equation will be exactly satisfied for a situation in which $ is uniform and 
(pU) varies linearly. With this additional term, the quantity Jl^ in equa- 
tion (57) is redefined as 

hi = Bi |(Ji-J 1+ i) - [(pU)i - (pU) i+1 ] 

+ C i |(Jj - J i _ 1 ) - [(pU) i - (pU) 1 _ 1 ] <f i _ 1 | (58) 

Note that the spline contribution Jl^ to the total flux is based on the dif- 
ferences in the J values at the adjacent faces of a control volume. Thus, a 
difference in J indicates the presence of a source and/or multidimensionality 
(a change of flux in one direction is felt as a source term in another direc- 
tion) . 


3.1.2 Discretization Equation for ft 

The governing equation (39) can be integrated over the control volume around 
the point i to yield 


J i+1 - J i = s Axi 

(59) 

The fluxes and Ji+i are replaced by the corresponding 
the linear flux-spline scheme, e.g., equation (57). The 
tion equation for 4> is 

expressions from 
resulting discretiza- 

ap = a E + a w 4>i_i + S c Ax* + S 

(60) 

where 


*E - H i+I exr < " P 1U ) 

(61) 

*w = "i “» < p i-i> 

(62) 

ap = a E + a w — Sp Axj 

(63) 

A /S A 

S = (Jli-Jl i+ i) 

(64) 
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In this last derivation, a linearized form of the source term S has been as- 
sumed, that is 


S = S c + Sp 

with 


( 65 ) 


Sp < 0 

The influence coefficients ag and ay appearing in the discretization equa- 
tion resulting from the linear flux-spline formulation are identical to those 
obtained from the exponential scheme. The contribution of the flux-spline is 
contained in the term S. 


Extension to Multi-Dimensional Situations 


In this section the extension of the linear flux-spline scheme to multidimen- 
sions is presented via the two-dimensional situation. The details for three- 
dimensional flows are very similar. 


For a two-dimensional situation, the conservation equation for the control vol- 
ume around the point (i,j) (Figure 11) is: 


+ (Jyi.j+l - J yij> Ax = S Ax Ay 


( 66 ) 



TE84-1681A 


Figure 11. Control volume in two dimensions. 
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The fluxes Jx and Jy are calculated using the one-dimensional flux-spline for- 
mulation in x and y direction, respectively. Thus, Jx and Jy are given by 


JX 1 ,J 

= (DXj 

*1-1, j E *1 *l,j > 

+ 


(67) 

Jy i,j 

= (Dyj 

♦i.J-1 ' Ey J *U ) 

+ 

31y l,J 

(68) 

final 

expression for $ is: 




•p *1, 

J “ a E 

*i+l,j + *W ♦i-l.J 

+ 

a N *i,j+l + *S *i,j-l 



+ s c 

Ax Ay + S 



(69) 


where 

a E = Hx i+1 exp (-Pxl +1 ) Ay (69a) 

a w = Hxi exp (Px^) Ay (69b) 

a N = Hyj+1 exp (-Pyj+i) Ax (69c) 

®S = H yj exp (PyJ_i) Ax (69d) 

a P = a E + a w + a N + a s - S P Ax Ay (69e) 

S = (Jlxi - Jlx 1+1 ) Ay (69f ) 

+ (Jiyj - Jlyj+i) Ax 


In equation (69), the descriptors x and y have been added to various quantities 
to denote the appropriate directions. These quantities can easily be obtained 
by substituting the appropriate directional parameters in the expressions de- 
rived earlier for one-dimensional flows. 


3.1.4 Solution Procedure 

A two-dimensional situation is governed by three field variables, $, Jx, and 
Jy. These variables are governed by the following three sets of equations: 

1. conservation equation for 4* 

2. spline continuity condition in the x-direction 

3. spline continuity condition in the y-direction 

Note that the spline-continuity conditions are one-dimensional in nature de- 
spite the multidimensional character of the whole problem. 

The relationship between the dependent variable ^ and its fluxes Jx and Jy 
is very similar to that encountered between pressure and velocity components 
in fluid flow calculations. The spline-continuity conditions give equations 
for the evaluation of flux fields. These equations involve $. The 4>-field 
should be such that the resulting fluxes satisfy the conservation equation for 
4>. Thus, the fluxes are analogous to the velocity components and the role 
of 4> is similar to that of pressure. Further, it should be noted that the 
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derivation of the ^-equation from the conservation of $ and the flux equa- 
tions is parallel to that of pressure in SIMPLER algorithm by combining the 
continuity and momentum equations. Based on these similarities the following 
procedure can be used for a coupled solution of Jx, Jy, and <J>. 

1. The flux terms Jlx and Jly are set equal to zero and the <J>-equation is 
solved. This solution will be identical to that from the (lower-order) 
exponential scheme. 

2. The computed $-field is used to obtain the fluxes Jx and Jy. 

3. These flux fields are used to compute 3lx and Jly. 

4. The (^-equation is solved with the spline contribution to the source 
term. 

5. New Jx and Jy fields^are obtained from the new $-field and the currently 
available values of Jlx and Jly using equations such as (57). 

6. Steps 3 through 5 are repeated until convergence is achieved. 

Note that a field solution for the fluxes is not performed. Instead, the 
fluxes are updated in an iterative manner (Jacobi update) • 

3.1.5 Fluid Flow Calculations 

The derivation previously presented is valid for the momentum equations also. 
Apart from the appearance of an additional source term due to the linear flux 
assumption, the discretization equation is similar to that from a lower— order 
formulation. 

The solution of the fluid flow equations involves a means of coupling the con- 
tinuity and momentum equations. These steps are identical to those used for a 
lower-order formulation and hence are not presented here. 

3.1.6 Notes on Computer Implementation 

1. The expressions for coefficients in the linear flux spline formulation in- 
volve the following two functions of the Peclet number: 

A(P) = 

e -1 

G(P) = £ CP-1) i 1 

If these expressions are used as such, the influence coefficients for 
linear flux-spline scheme are identical to those from the exponential 
scheme. To avoid the evaluation of computationally expensive exponentials, 
these expressions have been approximated by the following algebraic rela- 
tions: 

A(P) = max [0, (1 - 0.1 I P I) 5 ] + max [-P, 0] (70) 

G(P) = Q (P)/A (P) ( 71 ) 
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where 


Q (P) = 0.5 [1 + SGN (P)] - SGN(P) Q (- P ) (72) 

and 

Q ( _ P ) _ Q,5 - (O.U p ) 2 - (Q,?5454 P ) 6 

1 + P /6 + (0.25454 P )6 p 

With these new relations, the flux-spline scheme reduces to the Power-law 
differencing scheme under the assumption of constant flux within a control 
volume. 


2. The derivation of the flux-spline scheme is based on the assumption of a 
uniform diffusion coefficient within a control volume. The diffusion coef- 
ficients are stored at the main grid points (scalar locations) and are as- 
sumed to be constant over the control-volume around the grid point. Due 
to the staggered location of the velocity components, an estimate of the 
diffusivity over a velocity control volume is required. In addition, an 
interpolation is needed to evaluate the mass fluxes at the faces of the 
velocity control volumes. The practices used in the present implementation 
are as follows. Consider the U-control volume shown in Figure 12. 


The diffusion coefficient for is taken to be: 

+ 


r u 


r i..t tA Vi 


The mass flux at the west face is 
(pU) w = [(pU)i j + (pU) i _ 1>J ]/2 


(74) 


(75) 



TE88-4503 

Figure 12. Definition of a u-control volume. 
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and the mass flux at the south face is 


(pV) s - 


to i-i l Wu 


ta i-l + ta i 


( 76 ) 


3.2 CUBIC FLUX-SPLINE SCHEME 

The cubic flux-spline scheme is an improvement over the linear flux-spline 
scheme presented in Section 3.1. Similar to the linear flux-spline scheme, 
the details of the cubic flux-spline scheme are first presented for one-dimen- 
sional convection-diffusion situation. This is followed by the details related 
to its extension to multidimensions. 

3.2.1 One Dimensional Convectio n Diffusion 

In cubic flux-spline scheme, the total (convection + diffusion) flux (J) and 
mass flux (pU) are assumed to vary according to a cubic profile within a con- 
trol volume. For the control volume around point i, these are given by. 


J = 3 i + 


J ur J i 

Ax. 


X X 


^ + [ J ± ( J i+1 ^AXj ^AXj^ ^ 

3 2 

+ ^ J i +J i+1^ Ax i -2 ^ J i+l -J i^ ^Ax^ “ ^Ax^ ^ 


(77) 


pU = (pU) i + (pU)j x + |3[(pU) i+1 -(pU) i ] 

2 i 

- [2 ( P U)[ + (pU)^] AXjj (^) + {[(PU) ± + <P u > i+ i ] 


2 [(pU). , - (PU) 1 1 <5T) 


(78) 


where J*4 and (pU)'i are the derivatives of J and (pU), respectively, at the 
control-volume face. For ease of presentation, these equations are rewritten 

as: 

J = PO + Pl x + P2 x2 + P3 x3 (79 

(pU) = q 0 + Ql x + Q2 x2 + «l3 x3 (8 ° 

Substituting equations (79) and (80) in the definition of total flux, J, 


J 


pu<(> - r 


d$ 

dx 


(81) 
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and integrating, with uniform diffusivity I, the following expression for 
$ results: 


4> = A - exp (Jp(x)dx) J (C 0 + Cjx + C 2 x 2 
+ C 3 x 3 ) exp (- ,fp(x) dx) dx 


( 82 ) 


where 


P(x) = D 0 + D 2 x + D 2 x 2 + D 3 x 3 

Cf = Pl /r, ) 

D i = <Ii/r, i = 0, 1,2,3 ( 


(83) 


(84) 


The constant A in equation (82) is evaluated using the condition for d> at 
the main grid point 


♦ (x = Axj) = 4» ± 


(85) 


The expression for total flux J is obtained by imposing the continuity-of-6 
condition at a control-volume interface 

*1 = *t-l ( 86 ) 

After considerable algebra, the final expression for J can be written as (Ref 
4) : 


J i = D i (4>i-l 8i-l - ♦igl 1 ) + J1 + J3 
where 

Jl = Bj (Jj - J( + i) + Cj (Jj - Ji_^) 

23 - [BZj ( ^ h - to' + C2j toj+j 

+ B3 i t-toj + < J i + 2 i + i>J to- 


+ C3, t- 




toi 


^ - (J i + K-i 


„ , ii=i . ^ _L, 

1 r i-i n i-i r i W 


-1 


(87) 


(88) 


(89) 


(90) 


(Ax~) D t 12, 
1 = r 4 Ax 4 Ilj 


(91) 
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( 92 ) 


c (taj-i) D t I2 1 _ 1 g i _ 1 

1 r i-l **1-1 I1 i-1 


Ax 

B2 = — ~ B - EP 


ta i-i 

C2 i * 7TT c i - “i 


Ax 


i-1 


B3 i = EP t - 


(Axp D t I4 i 
(Axp 2 r i U i 


(Ax^) D 4 14 


C3 i = 2EM i - 


i-1 


/Ax 


i-1 


< to i-!> r i-l n i-l 


l Ax 


i-1/ 


(Ax ) D 13. 

EP i - ' r 1 ' Ax 1 ' 

EM <AXj!i) 2 P A 13^ 
1= r i-l *Vl n i-l 


il = s io dn 


12 = ii / n io dn 

o 


13 = ii J“ n io dn 

o 


14 = ii S n io dn 
o 


Io = exp [-(an + bn 2 + cn 2 + dn 4 )] 
g = exp (a+b+c+d) 


(93) 

(94) 

(95) 

(96) 

(97) 

(98) 

(99) 

( 100 ) 
( 101 ) 

( 102 ) 

(103) 

(104) 
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(105) 


a l " 


(pU^ Ax. i 

i\ 


, Ax. Ax. 

b i ' ( " V >1 2 f f 

I <P«) i+1 - <P»), , 

C 1 -| 3 ZT - 2 ‘<P u >i ♦ <p°>i + i J} 


(Ax.) Ax 

3 Ax, rT 


d i - 


j[(pU)j + (pU)^] - 2 [ 


(pU) i+l~ (pU) 

Ax, 




(Ax~) Ax i 


4(Ax i )" T i 


a i-i = (pu) i-i - (2b i-i + 3c i-i + 4d i-i> 


i-l 


[(pU) i - (pU)] i _ 1 Ax i ^ 1 Ax^ 

= Ax, , 2 I\ , ~ lm 


i-l 


5 c 


i-l 


i-l 


-1.5ci_i (1 - fi_i) 


(1+ft .) 3 
-2d — _i=l 


Ax 


c i_i = - (pu)t . J 


l=jl 


a=i 


i-l 6 (1 + f^x) r^! 
2di_i (1 - ft_ 2 ) 


(106) 


(107) 


(108) 


(109) 


(110) 


(111) 
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( 

[(pU). - (pU)^], 


d i-i = pp u >U + 

(pujp - 2 taj _ i 1 


S-1 , 

ta Yi 

(112) 

4 (1 + fj_ x ) 2 



f + = ^i-l 

i-1 .4 

to i-l 


(113) 


The expressions in other directions can be derived by axis transformation. 

The integrals in equations (99) to (102) are evaluated numerically. These de- 
tails can be found in Ref 4. 

The expressions presented previously involve the derivatives of the total flux 
and the mass flux. These derivatives are evaluated by solving the spline rela 
tionship between a cubic function and its derivative. The derivative of the 
total flux is obtained from the equation: 


i-1 


Ax 


+ 2 


i-1 


hh 


Ax 


i J 


J' 


i±l 


Ax, 


, J i+1 ~ J 1 . 3 J i I J l=l (114) 

- J O + ° 2 

^i to i-l 

The relationship for (p)’ is obtained by replacing J by (pU) in equation 
(114). 


3.2.2 Discretization Equation in Two Di mensions 

The discretization equation for 4> * 8 now obtained by substituting the flux 
expression, equation (87), in the conservation equation for 4>, equation (66). 
The general 4> discretization equation is: 

a P <|>i = a E 4>i+i + a w 4*1-1 + a N 4»j+i + a s 4>j-i + b (n 


where 

«E = DX i+1 j 1/gxi+ij , a W = DX i,J «*i-l,J 
a N = DY i,j+l i/^ij+l » a S = DY i,j ^i.J-l 
ap = a £ + a W + a N + a S “ Sp ^*1 ^ Y j + ® 
b = Sq Axj Ayj + b a( j 

B = max [al w ,0] + max [-al E ,0] + max [al s ,0] + max [-al N ,0] 


(116) 

(117) 

(118) 

(119) 

( 120 ) 
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b ad = jc i,j Ax* Ay} 

+ [(JIXij + J3Xj j) - (JlX i+1J + J3X i+ljj )] Ayj + 

+ tCJIYi j + J3Y lfj ) - (J1Y 1J+1 + J3Y iJ+1 )] Axi 

+ {max[-al w ,0] + max[al E ,0] +max [-al s ,0] + max[al N ,0] } (121) 

aly = CpU)ij Ayj - (a w>(i>j) -a E> (122) 

al E = (pU)i+l,j Ayj - (a w ,(i+l,j)-a E , (i,j)) (123) 

als = (pV)ij Axj - (a s ,(ij) -a N> (ij_i)) (124) 

ai N = ( P V) lfJ+1 Axi - (a s>(i>j+1) -a N> (iJ) ) (125 ) 


where all definitions given for the x direction, equations (87) through (114), 
should be applied for a line of constant j and analogous expressions should be 
used for the y direction. 

In the derivation of the flux balance over a control volume, it is generally 
assumed that the flux at the midpoint of a control— volume face prevails over 
the entire face. Further the value of the source at the grid point is assumed 
to prevail over the entire control volume. These approximations reduce the 
accuracy of the cubic flux-spline scheme. To compensate for this adverse ef- 
fect on accuracy, an additional correction term JCi,j (Ref 4) is added in the 
final discretization equation. This term is evaluated by ensuring that the 
conservation equation is exactly satisfied at the main grid point and involves 
the fluxes and their derivatives. 

3.2.3 Solution Procedure 

The main steps in the solution procedure for the cubic flux spline scheme ap- 
plied to convection-diffusion problems are as follows: 

1. guess a 4> field 

2. set all fluxes and flux derivatives equal to zero 

3. set all mass flow rate derivatives in the x and y directions equal to 
zero 

4. solve equation similar to (114) with the TDMA algorithm to obtain the 
pU)» field 

5. solve an equation analogous to (114) with the TDMA algorithm to obtain 
the pV)’ field 

6. calculate all flux coefficients, equations (90) through (113) for the x 
directions and analogous expressions for the y direction 

7. form all neighbor coefficients, equations (116) and (117), where all 
Integral functions are numerically evaluated by the use of a six point 
Gauss Quadrature 

8. calculate J1X and J3X, equations (88) and (89), by using the available 

Jx, and J'x fields. Proceed in an analogous way for the y direc- 
tion. Form the source term b and the ap term for the $ equation 
using equations (118) through (125) 

9. solve the $ equation (115) by the line-by-line TDMA algorithm 
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10. calculate new fluxes for the x direction from equation (87), where all 
neighboring Jx and J'x values are assumed known; proceed in an analogous 
way for the y direction 

11. solve the J'x equation (114) by the TDMA algorithm and similarly solve 
for J'y 

12. check convergence (If convergence has not been attained, one must repeat 
the process by going back to step 4. If the flow field is given and if 
the coefficients are constant, repeat the process going back to step 8. 

If the coefficients are constant, the process must be repeated from step 8 to 
avoid recalculating the integral functions. 

3.3 CONTROLLED NUMERICAL DIFFUSION WITH INTERNAL FEEDBACK (CONDIF) 

CONDIF is a modified central differencing scheme (CDS) that eliminates the 
over— and undershoots associated with CDS but retains its accuracy even at high 
Peclet numbers. CONDIF modifies the CDS by introducing a controlled amount of 
numerical diffusion based on the local gradients. For most problems, the 
numerical diffusion can be adjusted to be negligibly low. 

Since the CDS forms the basis of CONDIF, the development begins with a discus- 
sion of CDS. This is followed by the essential features of CONDIF. Like the 
derivation of other schemes, the analysis begins with one— dimensional transport 
of a scalar. 

3.3.1 Central Difference Scheme 

The one— dimensional convection-diffusion equation is 

= S (126) 
dx 


where 

J = po* - r g (127) 

In CDS, the values at a control-volume face are obtained using a piecewise 
linear profile for $. Thus, the values of $ and (d$/dx) at the west face 
of the control volume around point i (see Figure 13) are 


Ax~ Ax 

& m * A + — — * — * A 

T W * 6x 1 *i-l Sx l 9 


d$ 


!i±i=i 

6x. 


w i 

The discretized conservation equation is given by 


J e - J w = S Ax 


(128) 


(129) 


(130) 
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Figure 13. A one-dimensional situation. 
Substituting for the fluxes, the final ^-equation can be written as: 


where 


a„ = 


. * a E 4>i+l + a W 4>i-l + b 

(131) 

F e P°e 


6x, . ~ _ 

i+1 2 

(132a) 

r pu 

w K w 

6x^ + 2 

(132b) 

a E + a W 

(132c) 

S Ax 

(132d) 


For P (=pu5x/T) > 2, the coefficients in CDS become negative resulting in 
over- and undershoots in the numerical solution. CONDIF is a modified central 
differencing scheme that retains the essential features of CDS but ensures 
positivity of the coefficients. 

Consider the coefficients ag and a^; these may be expressed as 

a E = a E ~ A e (133) 

and 

a W = a W ~ A w (134) 


36 


where 


a E = r e /(6x) i+1 + [ KpU) e | - (pU) e ]/4 

a W = r w /6x)i + [l(pU) w | + (pU) w ]/4 

A e -[|(pn) e |+ (pU) e J/4 

A w = [ | <pU) w | - (pU) w ]/4 

Note that a£ and ay, A e , and A w are all unconditionally positive. 

Equation (131) may now be rearranged as: 

( a E + a w) 4»i = a E 4>i+l + a W $i-l (139) 

+ A e ($i - 4>i+l) + A w (4>i - 4>i-l) + b 

In equation (139), the terms involving A e and A w need special consideration 
if these terms are treated implicitly, the original CDS is recovered and oscil- 
lations are encountered in the numerical solution. Due to the dominant con- 
tribution of these terms at high Peclet number, an explicit treatment (i.e., 
as a source term) is also not feasible. The CONDIF scheme restructures this 
equation based on the local gradients of the dependent variable. 

3.3.2 Alternate Representation of CDS 

In equation (139), the term involving A w may be expressed as 

A e (4>i - $ i+ i) = A e r (4>i_i - 4>i); $i ¥ 4>i-i (140) 


(135) 

(136) 

(137) 

(138) 


where 

R = (4>i - 4>i+i)/(<l>i-i - 4>i) 

Similarly, the term involving A w may be rearranged as 

A w (4>i - 4>i_i) = a w (4>i+i - ti)/R; 4»i+i ¥ 4»i 
The special cases when 4*i equals or *re discussed later. 

Equation (139) may now be written as: 

(a'E + a 'w) 4>i = a ’E $i+l + a 'w $i-l 

with 

a’g — a*E + Ay/R 
a'y — + A e R 


(141) 


(142) 


(143) 


(144) 

(145) 
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Equation (143), for all practical purposes, is identical to equation (139) and 
is another representation of CDS. Also all terms in Equations (143), except 
R, are unconditionally positive. 

3.« 3 . 3 Derivation of CONDIF Scheme 

An examination R in equation (141) shows that, with reference to Figure 13, it 
may be written as 

R ~ (8<|>/8x) e /(9$/8x) w (146) 

R is proportional to the ratio of the gradient of 4> at the control-volume 
faces e and w. Finding the conditions under which R may become negative is 
now possible. For a monotonically increasing and decreasing function, R would 
always be positive. It is only for a function going through an extremum (a 
maximum or a minimum) within the control-volume around P that R becomes 
negative. 

In most fluid dynamic calculations, the dependent variables exhibit a monotonic 
behavior at most of the grid points; it is comparatively at a few grid points 
that the variables go through an extremum. With this observation, the CDS is 
modified as follows. 


1. For all the grid points where the grid Peclet number is less than or 
equal to 2, the CDS is used in the form given by equations (131) and 
(132). 

2. For the case where the grid Peclet number exceeds 2, the CDS is used in 
the form given by equations (143) to (145); if the R-parameter is posi- 
tive, otherwise. 

3. The upwind scheme is used. 


The suggested modification of CDS ensures the positivity of coefficients and 
hence eliminates the over- and undershoots. Note that the influence coeffi- 
cients now become quasi-linear. This is because the R-parameter represents 
the ratio of the gradient of the dependent at the control-volume faces. 


The parameter, R, may assume very high values if sharp variations in the gradi- 
ent of the dependent variable exist. Conversely R may approach a value of zero 
if $ is locally constant ($j = or $i+i). These cases repre- 

sent either a singularity in the governing equation or a trivial solution ($ 

= constant). Therefore, a limit on the R-parameter is imposed so that 


_1 

R 

max 


i R. 


max 


(147) 


where Rj^* is an arbitrary number. It determines the amount by which the 
gradient of the variable may change from one grid point to another. The 
parameter Rj^ plays the role of introducing a controlled amount of numeri- 
cal diffusion in the numerical diffusion. 
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The modified scheme previously described is named CONDIF (Controlled fiumerical 
Diffusion with Internal feedback). The numerical diffusion comes from the 
limit imposed on the R-parameter; the feedback refers to the self-adjusting 
feature of the scheme based on the ratio of the gradients at the adjacent con- 
trol-volume faces. 
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IV. EVALUATION OF SELECTED SCHEMES 


The three discretization schemes were used to solve a variety of two-dimen- 
sional test problems. The test problems included transport of a scalar, lami- 
nar flows and turbulent flows. For each case the numerical results were com- 
pared with a reference solution. In the turbulent flow test cases, emphasis 
was placed on the differences between the various numerical results rather than 
assessing their accuracy against the experimental data. Such a practice has 
been followed primarily because the use of a turbulence model introduces addi- 
tional uncertainty in the numerical results and it is difficult to distinguish 
between the model-related errors and the numerical errors. Further, the dif- 
ferences between the experimentally and numerically realized boundary condi- 
tions make a thorough comparison between the experimental data and calculations 
more difficult. To circumvent this difficulty, fine-grid numerical solutions 
have been taken as reference solutions for some test cases. 

The test problems used are 

scalar transport (2-D) 

1. recirculating flow with a prescribed heat source 

2. transport of a step discontinuity in a uniform flow at an angle 

laminar flows (2-D) 

1. flow over a backward-facing step 

2. flow and heat transfer in a shear-driven cavity 

turbulent flows (2-D) 

1. flow over a backward-facing step 

2. flow in a shear-driven cavity 

laminar flows (3-D) 

1. flow in a shear-driven cubic cavity 

turbulent flows (3-D) 

1. annular jet-induced flow in a duct 

2. row of jets in cross flow 

The results for these test cases are presented in the following sections. 

4.1 TWO-DIMENSIONAL TEST CASES 

4.1.1 Scalar Transport 

4. 1.1.1 Recirculating Flow with a Prescribed Heat Source (Revisit ed) 

This test case was described in subsection 2.3.3, where the performance of 
various schemes was presented in terms of an error parameter . The results for 
this problem are now analyzed in greater detail to assess the accuracy of the 
selected improved discretization schemes. Specifically, these results reveal 
the response of the discretization schemes to a variable source. 

The ^-distributions at the vertical midsection (x = 0) of the square domain, 
obtained using a uniform 22 x 22 grid, are compared with the exact solution in 
Figure 14. At P = 10, both flux-spline schemes almost reproduce the exact 
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Figure 14. ^-profiles at the midsection of the cavity; recirculating flow 

with a prescribed source. 
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solution. CONDIF results are significantly more accurate than those from the 
Power-law scheme, which produces a diffused profile, indicating the presence 
of false diffusion. 

The effect of false diffusion in the Power-law solution is much more evident 
at P = 100. The flux-spline schemes still give very accurate solution. CONDIF 
results, even though superior to those from the Power-law scheme, exhibit sub- 
stantial smearing, indicating the need of excessive numerical diffusion to 
stabilize the central-difference scheme at this high Peclet number and in the 
presence of a variable source. 

— S tep Discontinuity in a Uniform Flow at an Angle 

This problem is concerned with the transport of a step change in $ in a uni- 
form velocity field directed at angle to the x-axis as shown in Figure 15. 

With the diffusion coefficient r ■* 0, 4> is transported purely by convec- 
tion; consequently, the step change in $ at the upstream boundary is simply 
convected in the flow direction without any smearing. Even though there is no 
source present, this seemingly simple problem is a severe test for a convec- 
tion-diffusion scheme. The flow is directed at an angle to the grid lines, 
and the gradient in the cross-stream direction is infinite across the line AB 
in Figure 15. Further, there is no physical diffusion in this problem, and 
any smearing of the sharp profile would be a direct consequence of the false 
diffusion. 

For the given boundary conditions, the problem was solved using a uniform 11 x 
11 grid. Solutions were obtained at various flow angles, 9=0, 15, 30, and 
45 deg. The calculated profiles along the vertical centerline of the domain 
are presented in Figure 16. These results show that the profiles from the 



Vertical line for which 
_ the values of 0 are plotted 

* TE88-4505 

Figure 15. Transport of a step change in $ in a uniform velocity region. 
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Figure 16. Convective transport of a step change in profiles the 
vertical nidsection of the computational domain. 
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improved schemes are less smeared than those from the Power-law differencing 
scheme. Thus, the improved schemes introduce smaller, though finite, numerical 
diffusion. The linear flux-spline scheme reproduces the nodally exact solution 
® = 45 deg. For this problem, the flux-spline schemes also suffer from 
lack of boundedness. These oscillations are an indication of the dominance of 
the "flux" source terms at high Peclet numbers. The CONDIF results, on the 
other hand, do not show any oscillation. 

The over— and undershoot problem is inherent with most higher order discretiza- 
tion schemes. Various bounding strategies have been proposed to eliminate the 
unrealistic oscillations in the numerical solution. However, there is little 
experience with such methods. 

The oscillations associated with the improved schemes would diminish as the 
computational grid is refined. Further, in most practical problems, the pres- 
ence of physical diffusion will suppress these wiggles. With these points in 
mind, the presence of oscillations in the solution for this purely convective 
flow was not considered a serious shortcoming of the improved schemes. 

The next set of test problems involves the calculation of fluid flow. This 
requires an algorithm for the evaluation of the pressure field. For the re- 
sults reported here, the SIMPLE algorithm or its improved variant SIMPLER were 
used. These algorithms are based on a two— point differencing for pressure. 

The cubic flux-spline scheme is based on a cubic profile for the mass flux 
(pU) and the flux of the dependent variable. To preserve the accuracy char- 
acteristics of the cubic flux-spline scheme, a higher order representation of 
the pressure distribution is required. This calls for further research and 
development work, which was not undertaken in the present program. Thus, the 
fluid flow calculations presented next have been made using the linear flux- 
spline and CONDIF schemes only. Also, the linear flux-spline scheme has been 
simply referred to as the flux-spline scheme. 

4.1.2 Laminar Flows 

4. 1,2,1 Laminar Flow over a Backward-Facing Step 

The laminar flow over a backward facing step was computed using three convec- 
tion-diffusion schemes — Power-law, flux-spline, and CONDIF. The configuration 
for this test case is shown in Figure 17. Results were obtained for two values 



Figure 17. Flow over a backward facing step. 
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of Reynolds number, 100 and 715. The Reynolds number is based on the mean 
velocity in the inlet channel and the hydraulic diameter of this channel. The 
computed results are compared with the experimental data of Haas (Ref 14). 

The expansion ratio H/h is 2. 

The computational domain extended from the inlet plane to eight step heights 
downstream for Re = 100 and 32 step heights downstream for Re = 715. A fully 
developed velocity profile was prescribed at the inlet, and streamwise diffu- 
sion was neglected at the downstream boundary. The results were obtained using 
a 22 x 22 grid for Re = 100, whereas a 32 x 22 grid was used for Re = 715. 

The grid spacing was uniform in the transverse direction and increased geomet- 
rically with an expansion ratio of 1.05 in the axial direction. 

Figure 18 shows the computed and measured axial velocity profiles at various 
streamwise locations for Re = 100. At this low Reynolds number, there is very 
small false diffusion, and consequently, all three schemes give results that 
are in good agreement with the experimental data. The results from the flux- 
spline scheme show slightly better agreement with the experiments. For this 
case, C0NDIF was run with both Rmax = 1 and 10. Both values of Rmax 8 ave 
similar results, indicating, as expected, that false diffusion is negligible. 

The results at Re = 715 are shown in Figure 19. Notice that even though the 
use of improved schemes leads to better agreement between the experimental data 
and predictions at stations near the inlet, there is considerable discrepancy 
near the experimentally measured reattachment point (x = 13.3 h) . Similar be- 
havior has also been noticed by other investigators (e.g., Ref 15) and is at- 
tributed to the deviation of the flow from two-dimensionality at this Reynolds 
number. However, it should be mentioned that the reattachment length predicted 
by flux-spline (10.6 h) and C0NDIF (8.2 h) is longer than that from Power-law 
(7.2 h). This indicates a reduction in the false diffusion with the use of 
higher— order schemes. The longer reattachment length predicted by the flux- 
spline scheme compared with C0NDIF implies that there is less numerical dif- 
fusion in the former. 

4. 1.2. 2 Flow and Heat Transfer in a Driven Cavity 

The flow in a square cavity, with a moving wall (Figure 20), is a commonly used 
test problem for assessing the accuracy of recirculating flow calculations. 

Here in addition to the fluid flow, heat transfer calculations were also made. 
The temperature boundary conditions used were the moving wall at a temperature 
of unity and the remaining walls at a temperature of zero. The Prandtl number 
of the fluid was taken as unity. 

The problem was solved for a Reynolds number, based on the wall velocity and 
the cavity height, of 400 on a uniform 22 x 22 grid. The computed results are 
compared with a fine-grid (82 x 82) Power-law solution, which is labeled "REF- 
ERENCE" in the subsequent figures. 

Figures 21 and 22 show the distributions of the u and v velocities along the 
vertical and horizontal midplanes of the cavity, respectively. The flux-spline 
solution on a 22 x 22 compares very well with the REFERENCE solution. In fact, 
the two results are of comparable accuracy. The CONDIF results are also sub- 
stantially more accurate than those from the Power-law differencing scheme. 
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Figure 18. (continued) 







Figure 19. (continued) 




Figure 20. A driven cavity. 


The temperature distributions along the midsections of the cavity are shown in 
Figures 23 and 24. The f lux-spline results are very accurate. For this 
problem, the CONDIF results are slightly inferior to those from the Power- 
law scheme (PLDS) . This behavior is rather surprising since, for all other 
problems considered so far, CONDIF results were more accurate than PLDS re- 
suits. This point needs further investigation and was considered beyond the 
scope of the present work* 

The flow in this test case is strongly recirculating. In addition, the pres- 
sure gradients, which appear as sources in the momentum equations, are also 
significant. The accuracy of the flux-spline results demonstrates the capa- 
bility of the scheme to respond to flow skewness and presence of sources. Use 
of improved schemes, which do not account for the presence of sources, such as 
the skew upwind differencing scheme (Ref 7), for this problem will lead to re- 
sults that are only marginally better than those from the Power-law scheme 
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Figure 21. The u-velocity profile at the vertical midsection of the cavity. 



Figure 22. The v-velocity profile at the horizontal midsection of the cavity 











4.1.3 Turbulent Flows 

4 .1.3.1 Turbulent Flow over a Backwar d Facing Step 

The first turbulent flow test case considered is the flow over a backward fac- 
ing step designated as Stanford case 0421 (Ref 16)* 

The geometry of the test case is shown in Figure 25. The computational domain 
extended from the inflow boundary, shown dashed in Figure 25, to 20 step 
heights downstream. Initially, calculations were made using plug P rof “ es ** 
the inlet. These results were used for studying the effect of grid refinement 
as well as for a preliminary evaluation of the two improved differencing 
schemes. This was followed by calculations with inlet conditions corr « s P^"_ 
ing to the experiment of Kim (Ref 17). In the experiment, velocity and turbu 
lence data were not measured at the same locations. Velocity pro es are 
available at the step and at a distance four step heights upstream. Some tur- 
bulence data are also available upstream of the step. To prescribe the con 
Sons at ^e step, the following approach was taken. For the Power-law scheme 
the region upstream of the step was included in the computational domain A 
stepwise approximation to the experimentally measured veloc ty pro e 
scribed at the inlet, and the kinetic energy and the dissipation rates were 
calculated as follows: 

k = 0.0045 U 2 

e = 3 k 1>5 /H 

where U is the mean velocity in the upstream channel. The predicted profiies 
at the step were used as inlet conditions for the flux-spline and CONDIF 
schemes. 



H = 0.0762 m 
h * 0.03B1 m 
!_! = 0.1524 m 
L? = 2.3368 m 
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Figure 25. 


Geometry for turbulent flow over a backward facing step, 
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Table II. 

Turbulent flo w over a backward facing step 

calculated r eattachment lengths Tx p /hl * . 


Gtid 

Inlet 

profiles 

Power-law 

CONDIF 

Flux-SDline 

32 x 32 

Plug 

4.4 

4.2 

4.6 

40 x 40 

Plug 

5.0 

4.5 

5.3 

57 x 57 

Plug 

5.2 


5.3 

57 x 57 

Measured 

5.3 

— 

5.7 

*Measured value 

= 7 £ 1 





Results were obtained using three grid configurations of varying fineness. 

hC jf irSt grid was a unif °rm 32 x 32 grid selected without regard to the flow 
gradients. The second grid was a 40 x 40 nonuniform grid with more grid points 
near the walls and in the shear layer. Finally, computations were made on a 
nonuni form 57 x 57 grid to ensure grid independence. 


The overall accuracy of the computations can be judged by comparing the pre- 
dicted reattachment lengths with the measured value (7 + 1 h). The calculated 
reattachment lengths are presented in Table II. The flux-spline scheme shows 
an improvement over the Power-law scheme. Further, the flux-spline results 
reach an asymptotic value more quickly than those from the Power-law scheme 
which requires further grid refinement. 


T Jf, °£ thC C0NDIF scheme results in a smaller reattachment length compared 
with the Power-law scheme. This behavior is contrary to that shown by other 
l r ?! ed Rhemes including flux-spline. The cause for this behavior shown by 
CONDI F needs further investigation. Consequently, the fine grid calculations 
were not made using CONDIF. 


The agreement between the calculated and measured results improves slightly 
W Jf n measure d profiles are prescribed at the inlet. Even then, consider- 
able discrepancies remain between the two sets of results. Because the flux- 
spline results did not change as the grid was refined from 40 x 40 to 57 x 57 
* e p ^ asent flux-spline results can be considered grid-independent. Some fur- 
ther diagnostic runs were carried out by varying the inlet turbulence kinetic 
energy and dissipation length scale, but these changes had no influence on the 
reattachment length. 


The underprediction of the reattachment length for this problem has been noted 
in several previous studies (e.g., Ref 16, 18) and probably reflects a defi- 
ciency of the two-equation turbulence model. 
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The results for this test case, even though not conclusive, demonstrate the 
advantage of the improved discretization schemes in the evaluation of a physi- 
cal model since they have the potential of providing a grid-independent solu- 
tion without requiring an excessively fine grid. 

4. 1.3. 2 Turbulent Flow In a Shear-Driven Cavity 

This test case was selected to evaluate the performance of the flux-spline 
scheme for a highly recirculating flow. 

The configuration for this flow is identical to that shown in Figure 20* Cal- 
culations have been made at a Reynolds number of 10^. Results are presented 
for three uniformly spaced grids, 22 x 22, 42 x 42, and 62 x 62. No comparison 
has been made with the available experimental data. Instead, emphasis has been 
placed on the differences between the flux-spline and the Power-law results. 

Figures 26 and 27 show the velocity profiles along the vertical midsection of 
the cavity using the two differencing schemes on various grids. The flux- 
spline results show steeper gradients near the (stationary) lower wall. The 
flux spline results on the 42 x 42 grid are essentially grid-independent and 
do not change as the grid is further refined. This is in direct contrast to 
the Power-law results that are more sensitive to grid refinement. 

The effect of grid refinement is much more pronounced on the turbulence quanti- 
ties. Figures 28 and 29 show the computed turbulence kinetic energy and turbu- 
lent viscosity profiles along the vertical midsection of the cavity (x = 0.5), 
respectively. Both flux-spline and Power-law results on three grids have been 
included. Now the effects of grid refinement and differencing schemes are very 
clear. It is seen that the turbulence quantities are substantially under- 
predicted if the Power-law scheme is used, even on the finest grid. The flux- 
spline results approach the asymptotic limit with much fewer grid points. As 
the grid is refined from 22 x 22 to 62 x 62, the turbulent viscosity resulting 
from the Power-law scheme increases almost seven times. A similar change in 
the grid for the flux-spline scheme causes the turbulent viscosity to change 
by less than 10%. Similar trends are observed in the behavior of kinetic 
energy. 

The drastic changes in the levels of turbulent viscosity predicted by the 
Power-law scheme can be explained by examining the transport equations for k 
and e. In the central region, the shear is very small, and, therefore, the 
generation terms in the transport equations are negligible. However, the dis- 
sipation (sink) terms are finite. To adequately balance the nonzero sink 
terms, an accurate representation of the transport terms is needed. The 
numerical diffusion introduced by the use of the Power-law scheme will lead to 
a higher value of e since the sink term in its transport equation is propor- 
tional to the square of e. Since e also appears as the sink term in the 
k-equation, a low level of k will be predicted. Further, the sink term in e 
is proportional to 1/k; this causes overprediction of the levels of e. This 
double-edged effect results in very low values of turbulent viscosity. 
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Figure 26. The u-velocity profiles along the vertical midsection of the cavity. 
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Figure 27. The u-velocity profiles along the vertical midsection of the cavity. 
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Figure 28. Turbulence kinetic energy profiles along the vertical midsection 

of the cavity. 
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Figure 29. Turbulence viscosity profiles along the vertical midsection 

(for legends see Figure 28). 




4.2 SUMMARY OF THE TWO-DIMENSIONAL TEST CASES 

The performance of the three improved differencing schemes was evaluated by 
solving six test problems and comparing the results with available reference 
solutions . 

For the scalar transport equations both flux— spline schemes produced very accu- 
rate results. These schemes also showed the evidence of lack of boundedness 
in the purely convective flow. The CONDIF results were also superior to those 
from the Power— law scheme. By adding a sufficient amount of numerical damping 
in the regions of steep gradients, CONDIF eliminated the wiggles in numerical 
solution. To maintain its accuracy for fluid flow calculations, the cubic 
flux-spline scheme would require a higher order representation for pressure 
also. This work was not undertaken, and the cubic flux-spline scheme was not 
used for fluid flow calculations. 

The improved schemes showed improvement over the Power-law scheme for fluid 
flow calculations in both laminar flow test cases. The linear flux-spline 
scheme results were more accurate than those from CONDIF. A rather surprising 
finding was the inferior performance of the CONDIF scheme for heat transfer 
calculations in the shear-driven cavity. However, in an earlier problem in- 
volving the transport of a scalar in a recirculating flow, CONDIF results were 
more accurate than those from the Power-law scheme. This behavior of CONDIF 
requires further investigation. 

For the turbulent flow test cases, the flux-spline scheme was superior to the 
Power— law scheme. For the flow over the backward facing step, the use of 
CONDIF resulted in a reattachment length that was shorter than that predicted 
by the Power-law scheme. This behavior is contrary to that shown by other im- 
proved schemes for this problem. 

For the fluid flow calculations reported here, the flux-spline results were 
obtained using the SIMPLER algorithm (Ref 1) and the CONDIF results were ob- 
tained with the SIMPLE algorithm (Ref 1). For both schemes, the algebraic 
equations were solved using a line-by-line tridiagonal matrix algorithm (TDMA). 
No convergence difficulties were encountered with either scheme. 

4.3 SELECTION OF A SCHEME FOR THREE-DIMENSIONAL FLOWS 

For all the test problems considered so far, the linear-flux-spline scheme con- 
sistently produced more accurate results. The computational molecule for this 
scheme involves only seven points in three dimensions. Further, due to the 
one-dimensional nature of the spline-continuity conditions (even for a multi- 
dimensional problem), the extension of linear flux-spline scheme to three di- 
mensions is relatively straightforward. Due to these attributes, the linear 
flux-spline scheme was selected for incorporation into a computer program for 
three-dimensional flows. 
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A. 4 THREE-DIMENSIONAL TEST CASES 


4,4,1 Laminar Flow 

The flow situation under consideration is shown in Figure 30. Due to symmetry 
considerations, the computational domain extended only half the cavity width 
in the lateral (z) direction. 

Computations have been made at a Reynolds number, based on lid velocity and 
cavity depth, of 400. Results have been obtained using a uniformly spaced 22 
x 22 x 12 (x, y, z) grid. The present numerical results have been compared 
with the solution of Ku et al. (Ref 19) obtained using a pseudo-spectral method 
(25 x 25 x 13 modes). This solution has been designated as "REFERENCE" in the 
subsequent figures . 

Figures 31 and 32 show the velocity profiles of the u-component on the vertical 
centerline and v-component on the horizontal centerline of the symmetry plane 
z = 0.5, respectively. The flux-spline results are in better agreement with 
the reference solution than the Power-law results. In particular, the flux 
spline results exhibit sharper peaks than the Power-law solution due to less 
numerical diffusion in the former. 

The flow patterns in the y-z plane at x = 0.5 obtained using the two discreti- 
zation schemes are shown in Figure 33. The flux spline scheme predicts a more 
intense flow in the lower section of the cavity, especially in the vicinity of 
the eye of the vortex. 



Figure 30. Three-dimensional cubic cavity. 
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a. Flux-spline 


b. Power-law 


TE88-4512 

Figure 33. Shear driven cavity: flow pattern in the y-z plane at x = 0.5. 


— Turbulent Three-Dlroens tonal Jet-Induced Flow in a Duct 

4. 4. 2.1 Three-Dimensional Annular Jet-Induced Flow in a Duct 

This test case is one of the experiments being conducted tinder NASA HOST Ele- 
ment B for turbulence model evaluation. The test rig, shown in Figure 34, 
consists of five annular swirling or nonswirling streams issuing into a duct 
of rectangular cross section. There is provision for radial jets at specified 
axial locations. For the present test case, all streams are nonswirling and 
there are no primary jets. Due to symmetry of the problem, the computational 
domain includes only one quarter of an annular stream. In the cross section, 
the computational domain is bounded by symmetry lines in the z-direction and 
by a symmetry line and a wall in the y-direction. 

In the streamwise direction, the computational domain extended from the inlet 
plane to four duct widths downstream. Experiments indicated rather small 
changes in the velocity distributions beyond this location. 

Computations were made on two grids, 22 x 17 x 17 and 37 x 27 x 27. In the 
first grid, which will be referred to as "coarse,” the grid spacing was uniform 
in the cross section. In the refined grid, to be referred to as "fine,” a 
finer spacing was used within and near the jet. In both cases, the grid spac- 
ing in the x-direction was finer near the inlet. 

This test problem was solved using the Power-law and flux-spline schemes. Like 
the previous turbulent flow test cases, the emphasis has been placed on the 
differences between the results from these schemes, and very limited comparison 
has been made with the experimental data. This has been done primarily because 
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Figure 34. Geometry for an annular jet-induced flow in a duct. 


there is still considerable uncertainty about the adequacy of the turbulence 
model for such flows. Thus, the disagreement between the numerical and experi- 
mental results may be due to the numerical inaccuracy or the inadequacy of the 
turbulence model. However, for a given turbulence model, k-e model in this 
study, the flux-spline results are more accurate. 

Figure 35 shows the axial velocity profiles on the z = 0 plane (containing the 
centerline of the annular stream) at selected streamwise locations obtained 
using the coarse grid. For this calculation, plug profiles were prescribed at 
the inlet plane because the grid in the cross section was too coarse to allow 
for a distribution of various quantities. The use of the flux-spline scheme 
results in sharper peaks and a longer recirculation zone at the center. At all 
streamwise locations, the profiles resulting from the Power-law scheme are more 
smeared than those from the flux-spline scheme. This trend indicates the pres- 
ence of excessive numerical diffusion in the Power-law solution. 

Figure 36 shows the comparison of the predicted centerline velocity distribu- 
tion with the experimental data. The flux-spline results are in better agree- 
ment with the experimental data at locations near the inlet. However, further 
downstream the velocity is considerably underpredicted. 

The results for the fine grid are displayed in Figure 37. For this calcula- 
tion, the experimentally measured profiles of axial velocity and kinetic energy 
were prescribed at the inlet. The trends observed in the coarse-grid solutions 
are also noticed in this grid. The flux-spline scheme is able to preserve the 
peaks and produces profiles that are less smeared. Both schemes respond to 
grid refinement, especially near the inlet. To ensure grid-independence, cal- 
culations on a still finer grid would be required. 




Figure 35. Axial velocity profiles at z = 0 plane, coarse grid. 



1.2 



Figure 36. Variation of the centerline axial velocity, coarse grid 

(for legend see Figure 35). 


Due to the limitations of the resources, a further grid refinement was not 
undertaken. However, based on the experience with two-dimensional flows, the 
flux-spline results are not expected to change significantly as the grid is 
further refined. 

Figure 38 shows the computed centerline velocity variation on the fine grid. 

The flux-spline results are in good agreement with the experimental data. Both 
sets of calculations, however, now predict a stronger recirculation zone. 

4. 4. 2. 2 Row of Jets in Crossflow 

This three-dimensional turbulent flow test case, based on the experiments of 
Khan (Ref 20), involves a row of jets injected normal to the main flow in a 
duct of rectangular cross section. The physical situation for the particular 
case tinder consideration is shown in Figure 39. The jet diameter is 0.0254 
m. The pitch (jet centerline-to-centerline distance) is four jet diameters 
and the height of the test section is also four Jet diameters. 

The symmetry of the flow in the lateral (z) direction allows the computations 
to be confined between the jet centerline and the centerline between the jets. 
In the axial direction, the computational domain extends from 5 jet diameters 
upstream of the leading edge of the jet to 24 jet diameters downstream of the 
trailing edge of the jet. In the y-direction, the computational domain extends 
from the floor of the wind tunnel to the roof, or four jet diameters. 

At the inflow boundary (x/D = -5.5), plug profiles were assigned for all vari- 
ables. The streamwise diffusion was neglected at the outflow boundary. The 
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Figure 37. Axial velocity profiles at z = 0 plane, fine grid. 




Figure 38. Variation of centerline axial velocity, fine grid. 


circular Jet was approximated by an equivalent area rectangle with the longi- 
tudinal side equal to the jet diameter. The profiles at the exit of the jet 
were also taken as uniform. These assumptions may influence the degree of 
agreement between the experimental data and numerical predictions. Since the 
primary objective of these calculations is to stress the accuracy consideration 
without analyzing the consequences of the turbulence model being used, these 
simplified boundary conditions are considered adequate. 

Computations were made on two grids of varying fineness — 32 x 17 x 12 (coarse) 
and 37 x 22 x 17 (fine). Both grids were nonuniform with finer grid spacing 
near the jet. 

The results on the coarse grid are presented in Figure 40 in the form of the 
axial velocity profiles at four streamwise locations. These calculations show 
that both differencing schemes yield results that are significantly different 
from the experimental data. The numerical results, however, indicate that 
flux-spline calculations contain less numerical diffusion and hence are less 
smeared compared with those from the Power-law scheme. This is evident from 
the profiles at x/D = 4. The differences between the predictions and the ex- 
perimental data may be due to the coarseness of the mesh used and/or the inade- 
quacy of the turbulence model. To reduce the numerical errors, calculations 
were repeated on a finer grid. These results are shown in Figure 41. 

The use of the fine grid reduces the differences between the results from the 
two differencing schemes. The smeared Power-law profile indicates the presence 
of excessive numerical diffusion even on this grid. Although differences be- 
tween the calculated and experimental results are significant, further grid 
refinement may still be necessary to ensure grid independence of the present 


66 





Top vl#w 


Uoo 





Planes of 
symmetry 


Side view 



Figure 39. Flow geometry for "a row of jets in crossflow." 


calculations. Calculations for the same geometry using the bounded skew upwind 
differencing scheme BSUDS2 (Ref 21) indicate that even a 74 x 42 x 22 grid may 
not be sufficient to achieve grid independence, although the use of such a 
fine grid was not considered in this study. 

The effect of grid refinement is more pronounced on the flux-spline results, 
especially at x/D = 4. The Power-law scheme responds slowly to the change in 
grid. Such behavior is to be expected due to higher order of accuracy for the 
flux-spline scheme. Even on the fine grid, the Power-law solution is dominated 
by the false diffusion that may completely overwhelm the physical diffusion. 

Even though there are considerable differences between the experimental results 
and calculations presented here, the effects of grid refinement and change of 
differencing scheme are very similar to those seen in an earlier study (Ref 
22) for the same test case. 

This test case revealed some interesting issues related to the evaluation of 
the differencing schemes as well as the turbulence models. These primarily 
concern the specification of the inlet boundary conditions for the turbulence 
quantities. If a diffusive scheme, such as the Power-law scheme, is used for 
discretizing the equations and if the flow is recirculating, i.e., there is 
sufficient internal generation of turbulence, the results are insensitive to 
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Figure 44. Structure of the coefficient matrix. 


The matrix equation (148) is nominally linear because the elements of matrix A 
depend on the velocity components themselves, and also the source term b may 
be function of the dependent variables. In the present implementation, the 
nonlinearities are handled using the successive substitution (Picard) tech- 
nique, in which the coefficients and the source terms (including those arising 
from the flux-spline formulation) are calculated from the values of the depen- 
dent variables from the previous iteration. 

The direct solution of equation (148) was accomplished by the use of the Yale 
Sparse Matrix Package, YSMP (Ref 13). It solves the linear system A$ = b by 
a sparse matrix variation of the LU decomposition procedure, where L and U are 
the lower and upper triangular matrices. YSMP selects the elements of the main 
diagonal of the matrix as the pivots during the decomposition process. It is, 
therefore, necessary that the original coefficient matrix not contain any zero 
element on the main diagonal. As shown in Figure 44, the absence of pressure 
in the continuity equation leads to zeros along the main diagonal. This diffi- 
culty can be overcome in several ways, e.g., by rearranging the equations. In 
the present study, small nonzero elements are introduced along the main diago- 
nal in the continuity equation, and equation (148) is recast in terms of update 
vector A4> as follows: 
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A(A$)=b-A <fr* 


(149) 


where A<(>=<{> - 4>* and is the value from the previous iteration or ini- 
tial guess. Now, equation (149) is replaced by: 

A’ (A$) = b-A$* (15 ° 

where A f is the perturbed matrix that contains small nonzero elements along 
the main diagonal. Since A* is not used on the right side, the converged solu 
tion of equation (150) represents the solution of equation (148) with the con- 
tinuity equation properly satisfied. An added advantage of solving the equa- 
tions in terms of changes (A$) is that it reduces the round-off errors. 


Since the repeated numeric factorization of the coefficient matrix A is an ex 
pensive process, considerable execution time can be saved by not decomposing 
the coefficient matrix at each iteration. Instead, the decomposed matrix from 
a previous iteration is used. In effect, equation (150) is replaced by 

A*(A4>) = b-A$* <151) 

where A* is the old matrix decomposed at an earlier iteration, but A is the 
current coefficient matrix. In each iteration of this type, the right side is 
updated and the solution is obtained by forward and backward substitution using 
the previously computed decomposition. The cost of such an iteration is only 
10-15% of a full decomposition iteration. 

5.2 TREATMENT OF THE TURBULENCE QUANTITIES 

The turbulent flow calculations in this study are based on the k-c model, 
which requires the solution of two additional partial differential equations, 
after every iteration of the flow equations, to provide a new viscosity field. 

The two turbulence equations are highly coupled and nonlinear, due to their 
source terms, and difficult to solve. Due to these features, a direct coupled 
solution of these equations does not prove satisfactory and causes instabili- 
ties. Vanka (Ref 11) solves the k and c equations in a decoupled line-by- 
line manner by holding the other variable fixed. In the present work a 
slightly different approach was taken. Instead of solving the equations along 
a line and marching through the domain, the equations are solved sequentially 
over the entire computational domain using a line-by-line TDMA. The equations 
are solved repeatedly (typically 7 or 8 times) with updated source terms until 
K and c are sufficiently well converged. Such a procedure is expected to 
readily account for the ellipticity of the recirculating flows. 

5.3 PERFORMANCE OF THE ALGORITHM 

The coupled solution approach combined with the flux-spline scheme was applied 
to laminar flows in a square cavity and a planar sudden expansion and turbulent 
flow over a backward facing step. The details of these test cases are given 
in Table III. The number of iterations required for convergence and the execu- 
tion times for the coupled approach are compared with those for the sequential 
SIMPLER algorithm in Table IV. 

The performance of the coupled algorithm is nearly independent of the Reynolds 
number and the mesh aspect ratio. For the test problems considered, a coupled 
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Case No. 

Table III. 
Test cases. 

Reynolds 

Grid 

How 

number 

(uniform) 

Laminar flows 

i 

Cavity 

400 

22x22 

2 

Cavity 

1000 

22x22 

3 

Cavity 

400 

40x40 

4 

Cavity 

1000 

40x40 

5 

Expansion 

400 

22x12 

6 

Expansion 

400 

22x22 

7 

Expansion 

400 

42x22 

Turbulent flows 

8 

Backward facing step 

5.6xl0 5 

22x22 

9 

Backward facing step 

5.6xl0 5 

42x22 


Table IV. 

Number of iterations required and execution times. 




of iterations 

Execution 

times* 

Case No. 

SIMPLER 

Direct 

SIMPLER 

Direct 

1 

62 

17 

18 

6 

2 

84 

30 

24 

8 

3 



130 

21 

4 



140 

40 

5 

106 

47 

16 

5 

6 

122 

48 

35 

10 

7 



152 

42 

8 

800 

39 

400 

80 1 

9 



224 

58 1 

1 Convergence 

criteria 

for the turbulent 

flow cases are 

not same. 


* IBM 3084 cpu seconds 


solution of the continuity and momentum equations reduces the execution times 
by a factor of 3-5 compared with the SIMPLER algorithm. These findings are 
similar to those in earlier studies (Ref 10-12) where a lower-order convection- 
diffusion scheme was employed. 

The combination of a coupled solution approach and an improved discretization 
scheme offers a twofold advantage over the use of a sequential approach with a 
lower-order discretization scheme. For comparable accuracy, an improved scheme 
requires fewer grid points than a lower-order scheme. In addition, on a grid 
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the use of a coupled solution approach reduces the execution times compared 
with the sequential solution algorithms. Thus, when the two are combined, a 
significant reduction in execution times results. The advantage is expected 
to increase at high Peclet numbers and on finer grids. 

5.4 PLANE-BY-PLANE MARCHING PROCEDURE FOR 3-D FLOWS 

The coupled solution approach for fluid flow calculations presented previously 
in a two-dimensional context can easily be extended to three-dimensional flows. 
Such a procedure, however, will require a large amount of computer storage for 
LU factorization due to the increased number of nonzero elements in the coef- 
ficient matrix. These memory requirements are beyond the capacity of the cur- 
rently available computers. Hence, an alternate strategy needs to be devised 
for the extension of the coupled solution approach to three dimensions. In the 
present study, the equations are solved in a plane-by-plane manner. 

The predominant flow direction is selected as the marching direction. The 
equations in the planes normal to the marching direction are solved in a 
coupled manner using YSMP. There are two options available as to which vari- 
ables should be solved simultaneously. The three velocity components and pres- 
sure can be solved in a coupled manner or only the cross-stream (in-plane) 
velocities and pressure are solved implicitly and the axial velocity is solved 
separately. The first approach involves more nonzero elements in the coeffi- 
cient matrix and hence the cost of factorization is high. In the present im- 
plementation, the second approach has been followed; however, note that no com- 
parative study was conducted to evaluate the performance of these two alternate 
strategies. 

In the present plane-by-plane solution procedure, the treatment of in-plane 
velocities and pressure is identical to that in a two-dimensional problem. 

The additional features are related to the solution of the velocity component 
along the marching (axial) direction. It is necessary to ensure that the axial 
velocity component satisfies the axial momentum and continuity equations simul- 
taneously. This has been accomplished by following the procedure described in 
Ref 25. The essential steps in this method are as follows. 

The axial momentum equation for any cross-stream plane can be written as: 


a P u p 


E a nb ^nb + ^ ” 



(152) 


where nb refers to the neighbors in the cross section and Av is the volume 
of a control volume. 


In equation (152), b includes the physical source, the net transport in the 
axial direction and the axial pressure gradient given by the currently avail- 
able pressure field. (3p/3x) is the correction that must be applied to 
the calculated pressure gradient so that the axial velocity field at any plane 
yields the correct cross-sectional flow rate M: 

M = S p Up dA < 153 > 

The pressure gradient (3p/3x) should be such that both equations (152) and 
(153) are simultaneously satisfied. This is accomplished as follows. 
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For a fixed set of coefficients, there is a linear relationship between Up 
and (3p/3x). At a given cross-sectional plane, equation (152) can be writ- 
ten as: 


U 


P 




(154) 


The constants Up and dp can be evaluated by solving equation (152) with 
two different values of (9p/9x). With known values of these parameters, 
the required pressure gradient is calculated by enforcing equation (153). 
That is 


M=/ p U p dA = / pU p dA + J p dp dA 


(155) 


or 


© M - / p U d dA 

= Jpd dLA^ < 156 > 

P 

The final U-field is obtained by substituting the calculated value of (3p/3x) 
in equation (154). The pressure values at all locations downstream of the 
plane under consideration are also corrected using the calculated pressure 
gradient . 

The plane-by-plane solution procedure is supplemented by the solution of a 
three-dimensional pressure equation. This pressure field provides a means by 
which the elliptic effects are transmitted over the entire three-dimensional 
domain. In the present implementation, the pressure equation is identical to 
that in the SIMPLER algorithm. The pressure equation is solved after each 
sweep of the computational domain. 

5.5 EVALUATION OF THE PLANE -BY-PLANE SOLUTION TECHNIQUE 

The solution procedure described previously was applied to two model laminar 
flow problems. For these calculations, the Power-law differencing scheme was 
used. The flow situations considered are sketched in Figure 45. They are the 
following: 

1. flow in a three-dimensional sudden expansion 
2. flow in a rectangular box with a transverse jet 

These situations represent, in a simplified manner, the flow in a gas turbine 
combustor. 

The calculations presented here were started from zero initial guess for the 
velocity components and pressure. 

5.5.1 Three-Dimensional Sudden Expansion 


For the flow in a three-dimensional expansion, a 4:1 area ratio is considered. 
Due to symmetry considerations, computations have been restricted to a quarter 
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section of the duct. Two values of Reynolds number, defined as Re = Uin H/v, 
of 100 and 1000 are considered. Figures 46a and 46b show the convergence rate 
of the solution procedure for Reynolds numbers of 100 and 1000, respectively, 
on a uniform 22 x 12 x 12 grid. In these figures the absolute sum of residuals 
at all internal points for each equation are plotted against the number of 
iterations, which refer to the number of sweeps of the computational domain. 

5.5.2 Flow in a Box with Transverse Jets 

This flow was computed for a Reynolds number, Re = Vin H/v, of 500 using a 
uniform 22 x 12 x 12 grid. The convergence behavior of the solution procedure 
is shown in Figure 47. 

Both the flow situations considered here involve significant recirculation 
zones similar to that in a gas turbine combustor. On the grid used, the 
plane-by-plane procedure leads to convergence in 10 to 15 iterations. This 
indicates that the proposed procedure is robust. This solution technique was 
also used for solving the 3-D jet in cross flow turbulent test case described 
in section 4.4.2 on the coarse grid. Again, convergence was achieved in about 
20 iterations. In all these calculations the Power-law scheme was used. 

5.6 COST COMPARISON 

The plane-by-plane solution procedure is rapidly convergent and robust. A cost 
comparison with the iterative SIMPLER algorithm, however, revealed that the 
proposed procedure was not cost-effective, at least on the grids used. This 
high cost of solution is attributed to need to factorize the coefficient matrix 
at each cross-sectional plane, which is necessary to the rapidly changing na- 
ture of flow in the axial (marching) direction. Under specialized conditions, 
the factorized coefficient matrix at an upstream plane can be used for computa- 
tions at a particular plane. Such a practice, however, cannot be used for a 
general flow. 

5.7 NOTE ON THE USE OF FLUX-SPLINE SCHEME 

The objective in this program was to combine the plane-by-plane strategy with 
the flux-spline scheme. Initial exploratory rims indicated that for sample 
problems considered, such a combination did not lead to convergence except on 
very coarse grids. It should be mentioned that no such difficulty was encoun- 
tered in two-dimensional flows where all equations were treated in a fully 
coupled maimer. The only new feature in three-dimensional flows is the use of 
the plane-by-plane solution strategy. It is therefore likely that the decou- 
pling in the marching direction is responsible for the lack of convergence when 
the flux-spline discretization scheme is used. 
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Figure 46. Rate of convergence for 3-D sudden expansion 
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Rate of convergence for a jet in a rectangular box. 



VI . CONCLUDING REMARKS 


6.1 SUMMARY OF THE PRESENT WORK 

The objective of the present program was to develop accurate and efficient 
numerical schemes for predicting the flow field in gas turbine combustors. 

To improve the numerical accuracy various discretization schemes available in 
literature were evaluated on the basis of several criteria. Three schemes, 
linear flux-spline, cubic flux-spline, and CONDIF, were incorporated into a 
computer program for two-dimensional flows and were used to solve a series of 
test problems. The test problems included scalar transport, laminar flows, 
and turbulent flows. The numerical solutions were compared with available 
analytical solutions, experimental data, or fine grid numerical solutions. 

For all problems considered, the linear flux-spline scheme was consistently 
superior to other schemes and was selected for incorporation in a computer 
program for three-dimensional flows. 

The linear flux-spline scheme was used to solve several three-dimensional 
flows. For a given number of grid points, the flux-spline scheme produces re- 
sults that are far superior to those from the (lower-order) Power-law differ- 
encing scheme. Further, it has the potential of providing a grid independent 
solution without requiring an excessive number of grid points. 

To improve the computational efficiency of the overall solution procedure, a 
coupled solution approach for the fluid flow equations was adopted. In this 
approach, the discretized continuity and momentum equations are solved simul- 
taneously using a sparse matrix inversion technique. The coupled approach, 
when used in conjunction with the linear flux-spline scheme, reduced the execu- 
tion times by factors of 3 to 5 compared with the currently used sequential 
(SIMPLE-based) algorithms. 

The memory restrictions of the present computers do not allow a coupled solu- 
tion of all fluid flow equations in three-dimensional flows. To circumvent 
this difficulty, a plane-by-plane solution strategy was developed. In this 
procedure the cross stream velocities and pressure in a plane normal to the 
predominant flow direction are solved in a coupled manner and the streamwise 
velocity is solved separately. The plane-by-plane technique worked satisfac- 
torily when used with the Power-law discretization scheme. However, conver- 
gence could not be attained with the flux-spline scheme except on very coarse 
grids. This aspect needs further consideration. 

6.2 RECOMMENDATIONS FOR FURTHER WORK 

The linear flux-spline scheme used in this work for majority of calculations 
is susceptible to spatial oscillations for highly convective flows on coarse 
grids. Such a behavior may cause difficulties when the scheme is used to dis- 
cretize the transport equations of the "always positive" quantities such as 
the turbulence kinetic energy. It would be desirable to develop a bounded 
scheme that is free of the oscillations. 

The CONDIF scheme performed satisfactorily except for the heat transfer problem 
in the shear driven cavity and turbulent flow over a backward facing step. 

The causes for this anomaly should be identified and the scheme should be 
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tested for more problems Involving fluid flow. An attractive feature of CONDIF 
is its simplicity and ease of programming. 

For the test problems considered here, the cubic flux-spline scheme performed 
very well. It should be extended to fluid flow equations. 

The causes leading to the lack of convergence of the plane-by-plane solution 
strategy when used with the flux-spline scheme should be identified. In addi- 
tion, the use of other solution algorithms, such as the multigrid techniques, 
should be explored. 
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